Programming for Data Analysis

In R and Python

part 2 - python

In [1]:
# %%
import pandas as pd
import pylab as plt
import seaborn as sbn
# %%
In [2]:
ls data
info.txt                UPLC_Plasma_ExpDes.txt  UPLC_Plasma_Raw.txt
UPLC_Plasma_Clinic.txt  UPLC_Plasma_QC.txt
In [3]:
cat data/info.txt
 Dear all,

Please find attached the data for the lectures. It includes both raw and QCed UPLC plasma glycan data, information about exp design and several phenotypes (age,sex,case/control). Replicated standards have "stand" label and duplicated samples have "_D" label. For this cohort we have completely randomized design without blocking (we did not have any information about sex,age,case/control before randomization).

I hope the above is useful to you. Please let me know if you have any questions.

Kind regards,

Frano
In [108]:
UPLC_Plasma_Clinic = pd.read_csv("data/UPLC_Plasma_Clinic.txt",
                                 sep='\t',
                                 decimal=',',
                                 index_col='Sample')
UPLC_Plasma_Clinic.head()
Out[108]:
Sex Age CaseControl
Sample
ID_0772 male 58.412048 0
ID_0773 female 58.795345 0
ID_0956 female 58.954140 0
ID_0589 male 58.970570 0
ID_0782 female 59.069130 0
In [109]:
UPLC_Plasma_Clinic.replace({'CaseControl': 
                                {0: 'control', 1:'case'}}
                          ).head()
Out[109]:
Sex Age CaseControl
Sample
ID_0772 male 58.412048 control
ID_0773 female 58.795345 control
ID_0956 female 58.954140 control
ID_0589 male 58.970570 control
ID_0782 female 59.069130 control
In [110]:
fg = sbn.FacetGrid(data=UPLC_Plasma_Clinic, 
                   row='CaseControl', 
                   col='Sex')
fg.map(plt.hist, "Age")
Out[110]:
<seaborn.axisgrid.FacetGrid at 0x7fbd392f12b0>
In [111]:
with plt.xkcd():
    fg = sbn.FacetGrid(data=UPLC_Plasma_Clinic, 
                       row='CaseControl', 
                       col='Sex', size=3, aspect=2)
    fg.map(plt.hist, "Age", histtype='stepfilled')

bootstrapping

we would like to verify how anomalous are certain "structures" that we observe.

One way of doing so it to use our data to recreate similar looking random data and see how rare it is to observe a specific difference.

In [112]:
male_mean_age = UPLC_Plasma_Clinic.query("Sex=='male'")['Age'].mean()
female_mean_age = UPLC_Plasma_Clinic.query("Sex=='female'")['Age'].mean()
female_mean_age - male_mean_age
Out[112]:
-0.9518601281863823

how uncommon it would to see that difference if the distributions were the same?

In [113]:
UPLC_Plasma_Clinic['Age'].sample(n=5, replace=False)
Out[113]:
Sample
ID_0013    83.033539
ID_0014    62.272415
ID_1098    77.489388
ID_0305    68.832306
ID_0136    62.839149
Name: Age, dtype: float64
In [114]:
clone = UPLC_Plasma_Clinic.copy()
clone['Age'] = clone['Age'].sample(n=len(clone), replace=True).values

clone_male_mean_age = clone.query("Sex=='male'")['Age'].mean()
clone_female_mean_age = clone.query("Sex=='female'")['Age'].mean()
clone_female_mean_age - clone_male_mean_age
Out[114]:
-0.5532268359243915

We're missing something...what?

We should define this as a function, there is already enough clutter!

In [146]:
def sex_age_difference(df):
    """calculate the age differences between the genders"""
    male_mean_age = df.query("Sex=='male'")['Age'].mean()
    female_mean_age = df.query("Sex=='female'")['Age'].mean()
    return female_mean_age - male_mean_age

def reshuffle_age(df):
    """create a cloned dataframe with reshuffled ages"""
    clone = df.copy()
    clone['Age'] = clone['Age'].sample(n=len(clone), replace=True).values
    return clone

sex_age_difference(reshuffle_age(UPLC_Plasma_Clinic))
Out[146]:
-0.4945889744543308
In [130]:
replicas = plt.array([sex_age_difference(reshuffle_age(UPLC_Plasma_Clinic)) 
                      for i in range(20)])
replicas[:3]
Out[130]:
array([0.05087429, 0.06026722, 0.68936891])
In [131]:
with plt.xkcd():
    plt.hist(replicas, bins=20)
    plt.axvline(sex_age_difference(UPLC_Plasma_Clinic), color='r')
In [133]:
sum(replicas<sex_age_difference(UPLC_Plasma_Clinic))/len(replicas)
Out[133]:
0.008

to be able to replicate this result, we would have to "fix" the random numbers!

In [137]:
print(UPLC_Plasma_Clinic['Age'].sample().values)
print(UPLC_Plasma_Clinic['Age'].sample().values)
print(UPLC_Plasma_Clinic['Age'].sample().values)
[74.15194702]
[75.09103394]
[63.85763168]
In [139]:
print(UPLC_Plasma_Clinic['Age'].sample(random_state=1).values)
print(UPLC_Plasma_Clinic['Age'].sample(random_state=1).values)
print(UPLC_Plasma_Clinic['Age'].sample(random_state=1).values)
[61.88911819]
[61.88911819]
[61.88911819]

data transform

In [34]:
UPLC_Plasma_ExpDes = pd.read_csv("data/UPLC_Plasma_ExpDes.txt",
                                 sep='\t',
                                 decimal=',',
                                 index_col='Sample')
UPLC_Plasma_ExpDes.head()
Out[34]:
Plate Column Row
Sample
ID_0430 1 1 A
ID_1159 1 1 B
ID_0125 1 1 C
ID_1142 1 1 D
stand_38 1 1 E
In [35]:
len(UPLC_Plasma_Clinic), len(UPLC_Plasma_ExpDes)
Out[35]:
(1235, 1424)

wait...

let's re-read the details of the data:

Replicated standards have "stand" label and 
duplicated samples have "_D" label

so, if we ignore the replicated, we should have the same index, right?

In [36]:
UPLC_Plasma_ExpDes.drop(UPLC_Plasma_Clinic.index).head()
Out[36]:
Plate Column Row
Sample
stand_38 1 1 E
ID_0152_D 1 1 F
ID_0384_D 1 1 G
ID_0309_D 1 2 A
ID_0595_D 1 3 H
In [37]:
UPLC_Plasma_ExpDes.loc[UPLC_Plasma_Clinic.index].head()
Out[37]:
Plate Column Row
Sample
ID_0772 9 2 H
ID_0773 10 11 C
ID_0956 13 9 A
ID_0589 1 9 A
ID_0782 2 1 F
In [38]:
shared_index = UPLC_Plasma_ExpDes.index.isin(UPLC_Plasma_Clinic.index)
UPLC_Plasma_ExpDes[~shared_index].head()
Out[38]:
Plate Column Row
Sample
stand_38 1 1 E
ID_0152_D 1 1 F
ID_0384_D 1 1 G
ID_0309_D 1 2 A
ID_0595_D 1 3 H
In [39]:
is_stand = UPLC_Plasma_ExpDes.index.str.startswith('stand')
is_duplicated = UPLC_Plasma_ExpDes.index.str.endswith('_D')
UPLC_Plasma_ExpDes[is_duplicated].head()
Out[39]:
Plate Column Row
Sample
ID_0152_D 1 1 F
ID_0384_D 1 1 G
ID_0309_D 1 2 A
ID_0595_D 1 3 H
ID_0187_D 1 5 D
In [40]:
UPLC_Plasma_ExpDes[is_duplicated].index.str.replace('_D', '')
Out[40]:
Index(['ID_0152', 'ID_0384', 'ID_0309', 'ID_0595', 'ID_0187', 'ID_0267',
       'ID_0969', 'ID_0748', 'ID_0165', 'ID_1140',
       ...
       'ID_0778', 'ID_0470', 'ID_1166', 'ID_1025', 'ID_0963', 'ID_0412',
       'ID_1060', 'ID_0552', 'ID_0582', 'ID_0537'],
      dtype='object', name='Sample', length=111)

We need to start working with out transform to see what is going on.

There are two main families of transformations that one needs to understand:

  • grouping
  • joining
  • pivoting

grouping divide the dataset in sub dataset based on certain properties, apply an operation and merge the result together (such as calculating the average age for each Sex).

joining merge two different tables.

pivoting is used to transform something like a tidy table in a more squared table. The inverse operation is called melt.

Join

The joining can be applied in four ways:

  • inner join
  • outer join
  • left join
  • right join
In [43]:
UPLC_Plasma_Clinic.join(UPLC_Plasma_ExpDes).head()
Out[43]:
Sex Age CaseControl Plate Column Row
Sample
ID_0772 male 58.412048 0 9 2 H
ID_0773 female 58.795345 0 10 11 C
ID_0956 female 58.954140 0 13 9 A
ID_0589 male 58.970570 0 1 9 A
ID_0782 female 59.069130 0 2 1 F
In [42]:
pd.merge(UPLC_Plasma_Clinic, UPLC_Plasma_ExpDes,
         how='left',
         left_index=True,
         right_index=True,
        ).head()
Out[42]:
Sex Age CaseControl Plate Column Row
Sample
ID_0772 male 58.412048 0 9 2 H
ID_0773 female 58.795345 0 10 11 C
ID_0956 female 58.954140 0 13 9 A
ID_0589 male 58.970570 0 1 9 A
ID_0782 female 59.069130 0 2 1 F
In [45]:
joined = UPLC_Plasma_Clinic.join(UPLC_Plasma_ExpDes)

Grouping

Divide data in various subgroups, apply an operation to each one of the subgroups, merge the results together

In [46]:
joined.groupby(['Sex', 'Plate'])['CaseControl'].count()
Out[46]:
Sex     Plate
female  1        20
        2        23
        3        20
        4        14
        5        20
        6        23
        7        22
        8        18
        9        11
        10       16
        11       18
        12       16
        13       19
        14       16
        15       13
        16       19
male    1        58
        2        52
        3        58
        4        56
        5        54
        6        60
        7        51
        8        58
        9        65
        10       62
        11       59
        12       62
        13       59
        14       69
        15       60
        16       64
Name: CaseControl, dtype: int64
In [47]:
joined.groupby(['Sex', 'Plate'])['CaseControl'].count().unstack()
Out[47]:
Plate 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Sex
female 20 23 20 14 20 23 22 18 11 16 18 16 19 16 13 19
male 58 52 58 56 54 60 51 58 65 62 59 62 59 69 60 64

Pivoting e melting

they are the less intuitive operations.

Needed to convert long tables into wide tables and back.

Together with the joins, they are the main methods to manipulate tidy data into the shape most appropriate for analysis.

Let's see few examples to clarify

In [87]:
fake_data = [('Jane', '2016/01/01', 10),
             ('Jane', '2017/01/01', 11),
             ('Jane', '2018/01/01', 12),
             ('John', '2016/01/01', 8),
             #('John', '2017/01/01', 9), # this information is missing
             ('John', '2018/01/01', 10),
            ]
fake_data = pd.DataFrame(fake_data, columns=['name', 'date', 'value'])
fake_data
Out[87]:
name date value
0 Jane 2016/01/01 10
1 Jane 2017/01/01 11
2 Jane 2018/01/01 12
3 John 2016/01/01 8
4 John 2018/01/01 10
In [88]:
fake_data.pivot(index='name', columns='date', values='value')
Out[88]:
date 2016/01/01 2017/01/01 2018/01/01
name
Jane 10.0 11.0 12.0
John 8.0 NaN 10.0
In [89]:
fake_data.pivot_table(index='name', 
                      columns='date', 
                      values='value',
                      fill_value=0)
Out[89]:
date 2016/01/01 2017/01/01 2018/01/01
name
Jane 10 11 12
John 8 0 10
In [96]:
pivoted = fake_data.pivot_table(index='name', 
                                columns='date',
                                values='value',
                                fill_value=0)
pivoted.reset_index().melt(id_vars='name')
Out[96]:
name date value
0 Jane 2016/01/01 10
1 John 2016/01/01 8
2 Jane 2017/01/01 11
3 John 2017/01/01 0
4 Jane 2018/01/01 12
5 John 2018/01/01 10

In all honesty, melt is probably the least intuitive command of the whole pandas package...

In [22]:
fake_data = [('Jane', '2016/01/01', 10),
             ('Jane', '2017/01/01', 11),
             ('Jane', '2018/01/01', 12),
             ('John', '2016/01/01', 8),
             ('John', '2017/01/01', 11),
             ('John', '2017/01/01', 12),
             ('John', '2018/01/01', 10),
            ]
fake_data = pd.DataFrame(fake_data, columns=['name', 'date', 'value'])
fake_data
Out[22]:
name date value
0 Jane 2016/01/01 10
1 Jane 2017/01/01 11
2 Jane 2018/01/01 12
3 John 2016/01/01 8
4 John 2017/01/01 11
5 John 2017/01/01 12
6 John 2018/01/01 10

abbiamo valori multipli per alcuni giorni, bisogna decidere come aggregare questi valori (se usiamo pivot_table devono per forza essere numerici)

In [23]:
fake_data.pivot_table(index='name', 
                      columns='date', 
                      values='value', 
                      aggfunc=plt.mean)
Out[23]:
date 2016/01/01 2017/01/01 2018/01/01
name
Jane 10.0 11.0 12.0
John 8.0 11.5 10.0
In [24]:
fake_data.pivot_table(index='name', columns='date', values='value', aggfunc=max)
Out[24]:
date 2016/01/01 2017/01/01 2018/01/01
name
Jane 10 11 12
John 8 12 10

A simplified version of pivot and melt are the stack and unstack, that works on the dataframe indices and columns.

In [16]:
fake_data = pd.DataFrame(plt.randn(4, 4),
                         index=[['x', 'x', 'y', 'y'], [1, 2, 1, 2]],
                         columns=[['a', 'a', 'b', 'b'], ['c', 'd', 'c', 'd']],
                        )
fake_data
Out[16]:
a b
c d c d
x 1 1.165504 2.267955 -0.967565 0.014499
2 -0.451932 -1.830350 -0.521789 0.065337
y 1 -0.784047 0.893234 -0.380278 0.750762
2 1.763780 -0.673566 1.017364 -1.701770
In [17]:
fake_data.stack()
Out[17]:
a b
x 1 c 1.165504 -0.967565
d 2.267955 0.014499
2 c -0.451932 -0.521789
d -1.830350 0.065337
y 1 c -0.784047 -0.380278
d 0.893234 0.750762
2 c 1.763780 1.017364
d -0.673566 -1.701770
In [18]:
fake_data.unstack()
Out[18]:
a b
c d c d
1 2 1 2 1 2 1 2
x 1.165504 -0.451932 2.267955 -1.830350 -0.967565 -0.521789 0.014499 0.065337
y -0.784047 1.763780 0.893234 -0.673566 -0.380278 1.017364 0.750762 -1.701770

going back to our data...

In [50]:
pd.pivot_table(joined,
               index='Plate',
               columns='Sex',
               values='CaseControl',
               aggfunc=pd.Series.count,
               margins=True,
              )
Out[50]:
Sex female male All
Plate
1 20 58 78
2 23 52 75
3 20 58 78
4 14 56 70
5 20 54 74
6 23 60 83
7 22 51 73
8 18 58 76
9 11 65 76
10 16 62 78
11 18 59 77
12 16 62 78
13 19 59 78
14 16 69 85
15 13 60 73
16 19 64 83
All 288 947 1235
In [142]:
pd.pivot_table(joined,
               index='Plate',
               columns='Sex',
               values='CaseControl',
               aggfunc=plt.mean,
              )
Out[142]:
Sex female male
Plate
1 0.500000 0.293103
2 0.260870 0.326923
3 0.200000 0.379310
4 0.500000 0.375000
5 0.350000 0.259259
6 0.478261 0.333333
7 0.227273 0.254902
8 0.222222 0.413793
9 0.545455 0.276923
10 0.375000 0.370968
11 0.222222 0.322034
12 0.250000 0.370968
13 0.263158 0.271186
14 0.125000 0.347826
15 0.230769 0.250000
16 0.421053 0.234375
In [52]:
pd.pivot_table(joined,
               index='Plate',
               columns=['Sex', 'CaseControl'],
               values='Age',
               aggfunc=pd.Series.count,
               margins=True,
              ).astype(int)
Out[52]:
Sex female male All
CaseControl 0 1 0 1
Plate
1 10 10 41 17 78
2 17 6 35 17 75
3 16 4 36 22 78
4 7 7 35 21 70
5 13 7 40 14 74
6 12 11 40 20 83
7 17 5 38 13 73
8 14 4 34 24 76
9 5 6 47 18 76
10 10 6 39 23 78
11 14 4 40 19 77
12 12 4 39 23 78
13 14 5 43 16 78
14 14 2 45 24 85
15 10 3 45 15 73
16 11 8 49 15 83
All 196 92 646 301 1235
In [63]:
joined.groupby(['Sex', 'Plate'])['Age'].aggregate([plt.mean, plt.std])
Out[63]:
mean std
Sex Plate
female 1 70.958112 5.445270
2 68.704818 6.245644
3 69.742779 6.081070
4 70.910726 6.760206
5 68.620260 6.586808
6 69.472130 5.720965
7 69.391575 6.443942
8 66.330672 5.636552
9 71.004417 5.843082
10 71.523271 6.865312
11 66.122747 5.785487
12 67.710302 6.358500
13 67.116251 5.360202
14 68.465606 6.297453
15 66.832202 6.032745
16 68.785908 5.930863
male 1 69.743917 7.457167
2 69.917390 5.872961
3 69.397814 6.059217
4 69.705241 6.135054
5 70.550612 6.714617
6 70.432763 5.912694
7 69.821315 6.960955
8 70.628290 6.753675
9 69.613141 6.047489
10 69.814974 4.856323
11 69.184423 6.619136
12 69.838820 6.843786
13 69.355561 5.442088
14 69.145294 6.380279
15 70.684144 5.844680
16 68.941051 5.189806
In [150]:
(joined.groupby(['Plate', 'Sex'])['Age'].
    aggregate([plt.mean, plt.std]).unstack())
Out[150]:
mean std
Sex female male female male
Plate
1 70.958112 69.743917 5.445270 7.457167
2 68.704818 69.917390 6.245644 5.872961
3 69.742779 69.397814 6.081070 6.059217
4 70.910726 69.705241 6.760206 6.135054
5 68.620260 70.550612 6.586808 6.714617
6 69.472130 70.432763 5.720965 5.912694
7 69.391575 69.821315 6.443942 6.960955
8 66.330672 70.628290 5.636552 6.753675
9 71.004417 69.613141 5.843082 6.047489
10 71.523271 69.814974 6.865312 4.856323
11 66.122747 69.184423 5.785487 6.619136
12 67.710302 69.838820 6.358500 6.843786
13 67.116251 69.355561 5.360202 5.442088
14 68.465606 69.145294 6.297453 6.380279
15 66.832202 70.684144 6.032745 5.844680
16 68.785908 68.941051 5.930863 5.189806
In [68]:
grouped = joined.groupby(['Sex', 'Plate'])['Age']
aggregated = grouped.aggregate([plt.mean, plt.std]).unstack(level=0)
aggregated.head()
Out[68]:
mean std
Sex female male female male
Plate
1 70.958112 69.743917 5.445270 7.457167
2 68.704818 69.917390 6.245644 5.872961
3 69.742779 69.397814 6.081070 6.059217
4 70.910726 69.705241 6.760206 6.135054
5 68.620260 70.550612 6.586808 6.714617
In [85]:
aggregated['mean', 'female'] - aggregated['mean', 'male']
Out[85]:
Plate
1     1.214195
2    -1.212573
3     0.344965
4     1.205485
5    -1.930352
6    -0.960633
7    -0.429740
8    -4.297618
9     1.391276
10    1.708297
11   -3.061676
12   -2.128518
13   -2.239311
14   -0.679689
15   -3.851942
16   -0.155144
dtype: float64
In [75]:
aggregated['mean'].head()
Out[75]:
Sex female male
Plate
1 70.958112 69.743917
2 68.704818 69.917390
3 69.742779 69.397814
4 70.910726 69.705241
5 68.620260 70.550612
In [80]:
aggregated.xs('female', level='Sex', axis=1)
Out[80]:
mean std
Plate
1 70.958112 5.445270
2 68.704818 6.245644
3 69.742779 6.081070
4 70.910726 6.760206
5 68.620260 6.586808
6 69.472130 5.720965
7 69.391575 6.443942
8 66.330672 5.636552
9 71.004417 5.843082
10 71.523271 6.865312
11 66.122747 5.785487
12 67.710302 6.358500
13 67.116251 5.360202
14 68.465606 6.297453
15 66.832202 6.032745
16 68.785908 5.930863

Assignment

  • create a user on bitbucket
  • create a private repository
  • DO NOT INCLUDE THE DATA IN THE REPOSITORY!
  • include me and Lennart as partecipant (egiampieri)

  • check differencies on males and females for glycan levels

  • DOT NOT USE INFORMATIONS ABOUT CASES AND CONTROLS
In [ ]: